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During  the  lithiation  of  a  Si  anode  from  pure  Si  to  fully  lithiated  alloy,  the  volume  expands  four  times 
and  modulus  varies  by  several  tens  of  times.  Thus,  the  Li-ion  diffusion  and  the  stress  evolution  can  be 
strongly  coupled,  which  may  play  a  significant  role  in  determining  the  anode  performance.  In  this  work, 
we  present  a  theoretical  study  of  the  fully  coupled  diffusion  and  stresses  in  a  nonequilibrium  Li— Si  system 
by  taking  into  account  the  effects  of  composition-dependent  modulus,  finite  concentration,  and  boundary 
constraint.  The  Li-ion  diffusion  and  induced  stresses  in  a  bilayer  Cu-coated  Si  anode  at  the  nanometer  scale 
is  examined  to  show  these  important  effects.  The  transient  stress-assisted  diffusion  problem  is  solved 
numerically  by  a  finite  difference  method,  whilst  the  stress  field  is  obtained  analytically.  It  is  shown  that 
the  modulus  variation  with  composition  plays  a  mild  role  in  the  Li-ion  diffusion.  In  order  to  account  for  the 
finite  concentration  effect,  a  nonlinear  flux  equation  is  introduced  that  describes  the  Li-ion  diffusion  over 
the  full  range  of  concentration  from  dilute  to  near-saturation  state  in  a  unified,  symmetric  manner.  The 
finite  concentration  effect  is  significant,  especially  during  the  early  delithiation  process.  The  boundary 
constraint  effect  is  found  to  play  an  intriguing  role  in  the  chemical  diffusion.  The  bending  stress  results 
in  a  resisting  force  to  Li-ion  flow  preventing  effectively  the  Si  anode  from  full  lithiation.  The  constraint 
effect  is  significant  for  a  wide  range  of  Cu  thickness. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Volume  changes  when  a  second  phase  diffuses  into  or  out  of 
a  host  material.  If  the  volume  change  occurs  nonuniformly  or  is 
under  a  boundary  constraint  (e.g.,  a  film  being  bonded  to  a  sub¬ 
strate),  stresses  are  built  up  in  the  material.  Since  the  stresses  raise 
mechanical  deformation  energy,  one  would  intuitively  expect  the 
diffusion  process  that  is  undertaken  to  minimize  the  system  free 
energy  to  be  affected  in  exchange.  In  some  cases,  this  chemome- 
chanical  coupling  effect  can  be  significant.  The  chemically  induced 
stresses  can  be  so  high  as  to  cause  severe  morphological  changes, 
including  fracture.  This  effect  has  been  held  responsible  for  struc¬ 
tural  breakdown  of  active  electrode  materials  and  rapid  fading  of 
electrochemical  device  performances.  In  the  wake  of  searching  for 
new  battery  electrodes  capable  of  high  energy  density  and  high 
discharge  rate,  the  literature  has  recently  seen  a  surge  of  research 
interest  in  diffusion  induced  stresses  and  fractures  [1-39]. 
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The  theoretical  study  of  chemoelasticity  may  be  traced  back 
to  the  early  work  by  Prussin  [40].  The  follow-on  studies  have 
applied  the  theory  to  examine  the  diffusion  induced  stresses  and 
consequent  damage  processes  including  dislocations  and  fractures 
[41-45].  Most  of  these  studies  have  only  considered  the  influ¬ 
ence  of  diffusion  on  stresses.  A  few  others  have  included  the 
mutual  coupling  effects,  attempting  either  to  improve  the  accu¬ 
racy  of  an  analysis  or  to  interpret  interesting  phenomena,  such 
as  pressure  diffusion  where  the  diffusion  processes  respond  to  a 
stress  field.  Among  the  recent  studies  of  Li-ion  batteries,  stud¬ 
ies  in  Refs.  [7,15,18,31,32,34,36,37,39]  have  only  considered  the 
diffusion  induced  stresses.  Others  [1,2,6,8,19,20,24,30,35,38]  have 
considered  the  mutual  coupling  effects.  In  these  studies,  insights 
have  been  gained  into  the  fracture  of  active  electrode  particles 
and  the  delamination  of  active  electrode  films,  especially,  the  size 
effect  in  these  critical  processes  [15,18,30,31,34,36,37].  Atomistic 
modeling  and  first-principles  calculations  have  been  carried  out 
to  understand  the  underlying  physics  in  the  alloying/intercalation 
processes  in  electrodes  for  Li-ion  batteries  [9,1 0,1 6,22,28,29,33,35]. 
Based  on  such  a  bottom-up  approach,  it  was  suggested  that  the 
modulus  of  the  Li— Si  alloy  follows  the  rule  of  mixtures  [28]. 
This  effect  quantified  by  first  principles  calculations  has  been 
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considered  in  a  continuum  mechanics  study  by  a  couple  of  groups 
[18,35].  Besides  all  the  aforementioned  theoretical  studies,  exper¬ 
imental  measurements  of  stress  and  modulus  [25-27]  and  direct 
observations  of  lithiation  processes  [21]  have  been  reported. 
Though  scarce,  they  provided  crucial  insights  into  the  Li— Si  alloy¬ 
ing  process.  The  experimental  measurements  of  Young’s  modulus 
of  Lii2Si7  [46]  and  L^Sis  [47]  altogether  with  existing  data  points 
of  pure  Si  [48]  and  pure  Li  [49]  confirmed  the  rule  of  mixtures  being 
applicable  to  the  Li— Si  alloy  system  in  the  polycrystalline  form  (see 
Fig.  2).  However,  the  magnitudes  of  the  Young’s  moduli  of  these 
alloys  significantly  differ  from  the  theoretical  predictions  [28]. 

The  present  study  is  motivated  by  an  attempt  to  fabricate  low¬ 
dimensional  Si  anodes  and  improve  their  electronic  conductivity. 
Low-dimensional  Si  anodes  are  important  to  improve  the  capacity 
of  Li-ion  battery,  and  key  to  applications  such  as  electrical  vehi¬ 
cles  that  require  high  power  density  and  high  energy  density.  Free 
standing  Si  nanorods/wires/tubes  have  been  fabricated  for  anode 
usage  in  Li-ion  batteries  [3,13].  Although  their  performance  can  be 
significantly  improved  relative  to  that  of  bulk  and  thin-film  Si  mate¬ 
rials,  they  suffer  from  severe  morphological  changes  over  cycles, 
intensified  solid  electrolyte  interphase  buildup  in  the  presence  of 
larger  surface/interface  area,  and  rapid  capacity  fading  compared  to 
other  types  of  anodes.  To  address  some  of  these  issues,  Cui  et  al.  [  1 1  ] 
proposed  a  carbon  nanofiber-Si  core-shell  structure  that  may  hold 
the  one-dimensional  active  material  together  and  connect  it  to  the 
current  collector  upon  fragmentation  (as  long  as  the  carbon  nan¬ 
otube  is  connected).  An  alternative  is  to  partially  coat  a  Si  nanorod 
with  a  Cu  film,  or  other  material  combinations  [50].  The  Cu  film 
would  provide  structural  support  to  the  Si  nanorod.  Meanwhile,  it 
can  improve  the  electronic  conductivity  of  the  anode  structure  and 
hence  reduce  thermal  issues. 

The  present  study  aims  to  examine  the  Li-ion  diffusion  and 
stress  evolution  in  a  Si  nano-anode  partially  supported  by  a  Cu 
coating.  The  Cu  film  would  impose  partial  constraint  on  the  Si 
nanorod  deformation  during  the  lithiation  and  delithiation  pro¬ 
cesses.  Because  Cu  is  inactive  to  Li-ion  alloying,  it  would  partially 
insulate  the  Si  nanorod  from  Li-ion  contact  and  cause  uneven  Li- 
ion  diffusion  on  the  opposite  surfaces  of  the  nanorod.  The  uneven 
Li-ion  diffusion  would  result  in  the  nanorod  bending,  which  may 
or  may  not  be  favored  for  the  purpose  of  stress  relief  without 
compromising  charge/discharge  rate  or  capacity.  During  the  lithi- 
ation/delithiation  process,  the  modulus  of  the  nanorod  varies  by  a 
couple  of  tens  of  times  due  to  the  composition  change,  and  the  Li- 
ion  diffusion  experiences  the  full  range  of  concentration  variation, 
from  dilute  to  saturated  Li-ion  concentration.  The  present  study 
theoretically  addresses  these  effects  of  composition-dependent 
modulus,  finite  concentration  and  boundary  constraint  on  Li-ion 
diffusion  in  a  Si  nano-anode.  Based  on  the  analysis,  the  following 
results  are  worth  remarking:  (a)  a  fully  coupled  chemoelastic  anal¬ 
ysis  of  Li-ion  alloying  with  Si  showing  quantitatively  the  effects  of 
composition-dependent  modulus  in  the  alloying  process;  (b)  intro¬ 
duction  of  a  nonlinear  flux  equation  that  describes  Li-ion  diffusion 
over  the  entire  range  of  concentration  in  a  unified,  symmetric  man¬ 
ner;  (c)  demonstration  of  intriguing  boundary  constraint  effect  in 
a  bilayer  Cu-coated  Si  anode. 


2.  Problem  formulation 

Consider  a  bilayer  Cu-coated  Si  anode  as  shown  in  Fig.  1. 
The  two  materials  are  perfectly  bonded  at  the  interface.  The  Si 
layer  is  subjected  to  Li-ion  alloying.  Since  Cu  is  inactive  to  Li-ion 
alloying,  it  blocks  Li  ions  from  entering  the  Si  layer  through  the 
covered  surface.  It  also  acts  as  a  structural  support  to  the  Si  layer. 
Both  of  the  diffusion  and  deformation  fields  are  assumed  to  be 
one-dimensional  through  the  thickness.  The  Cartesian  coordinate 
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Fig.  1.  Schematic  showing  a  bilayer  Cu-coated  Si  anode  that  is  subjected  to  Li-ion 
alloying.  The  diffusion  and  stress  fields  are  both  assumed  to  be  one-dimensional 
along  the  x-axis.  The  composite  beam  may  elongate  along  the  z-axis  and  bend  about 
they-axis. 


system  is  established  with  the  x-axis  directing  along  the  thickness 
direction  and  the  origin  located  at  the  left  boundary,  as  shown  in 
Fig.  1.  There  is  no  constraint  for  deformation  along  the  z-axis  or  for 
bending  about  the  y-axis. 


2.2.  Equations  of  chemo elasticity 


Within  the  classical  Euler-Bernoulli  beam  theory,  a  plane  per¬ 
pendicular  to  the  beam  boundary  is  assumed  to  remain  planar  and 
perpendicular  to  the  beam  boundary  upon  deformation.  It  results 
in  a  linear  variation  of  strain  throughout  thickness: 


£z  =  £zo  —  kyx,  (1) 

where  ez{=  duz/dz)  is  the  total  strain  (the  only  nontrivial  z- 
component),  uz  is  the  z  component  of  displacement,  sz0  is  the  strain 
at  the  boundary  (x  =  0),  and  ky  is  the  curvature  of  the  beam.  By  sub¬ 
tracting  the  contribution  due  to  inserted  Li  ions  (if  any),  the  elastic 
strain  field  is  given  by 

_  /  £zo  -  kyx  ~  yc  forO  <  x  <  hinSi,  .  . 

6z  ~  \  £zo~kyX  for/i  <  x  <  (1  +  n)/iinCu, 

where  h  is  the  Si  layer  thickness,  n  is  the  thickness  ratio  of  the  Cu 
layer  to  the  Si  layer,  c  is  the  Li-ion  concentration,  and  y  is  the  linear 
chemical  expansion  coefficient  characterizing  dimensional  change 
due  to  Li-ion  insertion  in  Si.  The  stress  field  is  given  by 


j  E{sz o  -  kyx  -  yc)  for  0  <  x  <  h  in  Si, 

1  Ecuteo  -  kyx)  for/i  <  x  <  (1  +  n)/iinCu, 


(3) 


where  E  and  ECu  are  the  Young’s  moduli  of  the  Li— Si  alloy  and  the  Cu 
layers,  respectively.  For  the  sake  of  brevity,  the  subscript  notation, 
Li— Si,  is  omitted  for  the  modulus  of  Li— Si  alloy.  The  alloy  modulus 
E(c)  varies  with  c,  which  will  be  described  later.  The  copper  modulus 
ECu  is  a  constant.  The  equilibrium  conditions  are  expressed  by 


*(1  +n)h 


(jzdx  —  Fz , 


*(1  +n)h 


azxdx  =  My  -  i Fzh , 
2 


(4) 

(5) 


where  Fz  and  My  are  the  force  and  couple  moment  applied  at  the 
centroid  of  the  Si  layer,  respectively.  Both  of  Fz  and  My  are  set  to  be 
zero  in  later  simulations. 
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Substituting  Eq.  (3)  in  Eqs.  (4)  and  (5)  and  rearranging  result  in 


(E  +  ECun)eZQ 


_  1  7 

Ex+  -£Cu/i(rr  +n) 


ky  —  yEc  +  Fz , 


(6) 


—  1  n 

£x  +  -£Cu/i(rr  +  n) 


£z0 


1 


Ex2  +  EECuh2{n 3  +  3n2  +  3n) 


/<v 


—  y  Ecx  ATy  . 


(7) 


In  the  above  equations,  the  over-bar  denotes  an  average  of  the 
covered  function  over  the  Si  layer,  for  example, 


Ecxdx. 


(8) 


Given  Fz,  My  and  c(x),  and  given  the  material  constants  and  geo¬ 
metrical  parameters,  Eqs.  (6)  and  (7)  can  be  analytically  solved  to 
determine  sz0  and  ky ,  with  which  the  stress  and  strain  fields  in  the 
composite  beam  can  be  fully  calculated. 


2.2.  Equations  of  mechanochemical  diffusion 


phase  has  been  observed  during  the  lithiation  of  a  crystalline  Si 
nanowire  [3,53]  that  does  not  seem  to  be  explainable  with  this 
assumption.  Although  the  underlying  physics  can  be  different,  this 
speculation  of  similar  composition  dependence  of  Li-ion  mobility 
during  the  Si  lithiation  and  delithiation  is  rooted  on  our  recent  study 
of  H  diffusion  in  MgHx  with  the  aforementioned  exponential  con¬ 
centration  dependence  of  H  mobility  and  modified  Fields  first  law 
given  in  Eq.  ( 1 0).  It  was  shown  that  a  sharp  phase  boundary  of  Mg 
hydride  can  arise  during  the  dehydrogenation  process  (but  not  dur¬ 
ing  the  hydrogenation  process)  as  a  natural  result  of  the  exponential 
concentration  dependence  of  H  mobility  [54]. 

By  substituting  Eq.  (9)  in  Eq.  (10),  the  flux  equation  is  rewritten 
as 

J  =  („> 

Carrying  out  the  derivative  in  the  second  term  on  the  right-hand 
side  and  rearranging  yield 

J  =  -D*^-Mc(\-c)(yE-Ifaf)ky,  (12) 

where  D  is  the  effective  diffusion  constant,  defined  as 


The  Larche-Cahn  chemical  potential  [42]  is  applied  to  the  Li— Si 
system: 


D*  = 


ycsE 

1  +  _c) 


2  E  c  E  cc  2 

y  -  -^r Oz  +  -^07 


E 2 


2yE3 


MRT, 


(13) 


/x  =  jlcq  +  RT  In 


c 


cs-c 


E  c  ? 
^+2E^’ 


(9) 


where  R  is  the  universal  gas  constant,  T  is  the  temperature,  cs  is  the 
saturation  concentration  at  the  stoichiometric  limit  of  Li22Si5,  and 
E?c  =  dE/dc.  The  first  term  /x0  is  invariant  of  composition.  The  sec¬ 
ond  term  represents  the  contribution  of  configurational  entropy, 
where  the  finite  concentration  effect  is  taken  into  account.  The  third 
term  is  due  to  the  volume  change  coupled  with  applied  stress.  The 
last  term  is  due  to  the  composition  dependence  of  modulus.  Larche 
[44]  suggested  that  this  last  term,  which  is  a  quadratic  function  of 
stress,  would  be  small  compared  to  the  preceding  term,  which  is  a 
linear  function  of  stress,  because  the  stress  is  typically  small  (rela¬ 
tive  to  E).  Sethuraman  et  al.  [27]  estimated  the  last  term  to  be  about 
only  1%  in  magnitude  of  the  preceding  term  for  the  Li— Si  system. 
Later,  we  will  show  that  the  last  term  can  make  a  mild  difference  up 
to  15%  in  the  predicted  stresses.  Therefore,  this  last  term  is  kept  in 
our  formulation.  An  interpretation  of  the  discrepancy  will  be  given 
later  with  respect  to  the  composition  dependence  of  E. 

The  Li-ion  diffusion  flux  in  Si  is  assumed  to  follow  a  nonlinear 
law: 


J  =  — Mc(l  —  c) 


dfi 

dx 


(10) 


where  M  is  the  Li-ion  mobility,  and  c(=  c/cs)  is  the  normalized 
Li-ion  concentration.  Compared  to  the  classical  flux  equation  relat¬ 
ing  to  chemical  potential  gradient  [51],  the  extra  modifying  factor 
(1  -  c)  considers  the  situation  where  only  those  lithium  ions  with 
immediate  neighboring  interstitial  sites  vacant  can  migrate.  It  is 
also  assumed  that  the  probability  for  a  lithium  ion  finding  its  imme¬ 
diate  neighboring  interstitial  site  vacant  is  uniform.  When  c  is  small, 
compared  to  cs,  Eq.  (10)  is  reduced  to  the  common  flux  equation 
by  neglecting  c  in  (1  -  c).  When  c  is  close  to  cSl  i.e.,  when  c  «  1, 
Eq.  (10)  is  then  reduced  to  the  common  flux  equation  for  the  dif¬ 
fusion  of  dilute  vacancies  at  concentration  (cs  -  c).  Recently,  Yang 
et  al.  [52]  identified  through  a  combined  experimental  and  theo¬ 
retical  study  that  the  mobility  of  interstitial  H  atoms  in  Mg  and  the 
mobility  of  H-vacancies  in  MgH2  can  differ  by  three  orders  of  mag¬ 
nitude.  In  the  present  study  of  a  Li— Si  system,  it  is  assumed  that 
the  Li-ion  mobility  is  independent  of  composition.  This  commonly 
adopted  assumption  should  be  subjected  to  further  examination 
because  a  sharp  phase  boundary  between  pure  Si  and  its  lithiated 


where  £cc  =  (92£/9c2).  The  mass  conservation  equation  is  given  by 

*  =  -1.  (14) 

dt  dx  K  J 

We  impose  a  boundary  condition  relating  the  boundary  concen¬ 
tration  and  flux  linearly: 


J  =A(ceq-c)  atx  =  0, 
J  =  0  atx  =  h, 


(15) 

(16) 


where  A  is  a  non-negative  rate  constant,  and  ceq  is  a  target  (equilib¬ 
rium)  concentration.  If  c  is  below  ceq ,  lithium  ions  flow  into  the  Si 
layer.  In  the  case  when  c  is  above  ceq ,  lithium  ions  leave  the  anode. 
In  the  process  of  lithiation,  ceq  is  set  to  be  the  saturation  concen¬ 
tration  cs.  In  the  reversed  process  of  delithiation,  ceq  is  switched  to 
zero.  This  condition  embraces  a  rate  constant  that  can  be  adjusted 
to  impose  different  charging  and  discharging  rates. 


2.3.  Materials  properties 


Most  recently  Ratchford  et  al.  [46,47]  measured  the  Young’s 
moduli  of  polycrystalline  Lii2Si7  and  Li22Si5  to  be  52GPa  and 
35.4  GPa,  respectively.  Based  on  the  vast  literature  of  Si,  the  Young’s 
modulus  of  polycrystalline  Si  may  be  assumed  to  be  170  GPa  [48]. 
Lastly,  the  Young’s  modulus  of  polycrystalline  Li  was  reported  to 
be  8  GPa  in  an  earlier  paper  [49].  These  four  data  points  are  plot¬ 
ted  versus  Li  atomic  fraction  x/(x+y)  in  LixSiy  alloy  (where  x  andy 
temporarily  used  here  for  atomic  fraction  should  not  be  confused 
with  axes  x  andy),  as  shown  in  Fig.  2.  A  linear  curve  fitting  is  found 
to  approximate  well  these  data  points,  suggesting  the  applicability 
of  the  rule  of  mixtures  in  this  case.  This  behavior  is  similar  to  what 
Shenoy  et  al.  [28]  suggested  based  on  their  first-principles  calcula¬ 
tion.  However,  it  should  be  remarked  that  their  predicted  moduli  at 
high  Li  concentrations  are  significantly  greater  in  magnitude  than 
the  experimental  results.  The  following  formula  of  composition- 
dependent  Young’s  modulus  is  used  in  later  simulations: 


£  =  8  + 


162 

4.4c  +  1 


(GPa). 


(17) 


It  recovers  exactly  the  two  data  points  of  Si  [48]  and  Li  [49] 
but  approximately  those  of  Li— Si  alloys  [46,47].  More  experimental 
data  points  are  required  to  affirm  Eq.  (17). 
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Fig.  4.  Schematic  showing  a  finite  difference  discretization  of  the  x  and  time  axes. 


Fig.  2.  Variation  of  LixSiy  alloy  Young’s  modulus  with  Li  atomic  fraction  x/(x+y).  The 
dashed  line  indicates  a  linear  fitting  curve  to  the  data  points  (symbols). 


Based  on  the  volume  and  space  groups  of  a  few  intermedi¬ 
ate  crystalline  phases  of  LixSiy  alloy  given  in  Table  1  in  Ref.  by 
Shenoy  et  al.  [28],  we  plotted  the  variation  of  alloy  volume  per 
Si  atom  as  a  function  of  composition  ratio  x/y,  as  shown  in  Fig.  3. 
It  shows  that  the  alloy  volume  is  nearly  a  linear  function  of  x/y. 
In  the  present  study,  the  Li-ion  concentration  c  is  defined  as  the 
amount  of  Li-ion  per  unit  initial  volume.  Thus,  c  is  proportional  to 
composition  ratio  x/y.  Based  on  this  graph,  the  chemical  expan¬ 
sion  coefficient,  y  can  be  reasonably  assumed  to  be  a  constant,  and 
is  estimated  to  be  1.64  x  10-6m3moL1.  The  saturation  concen¬ 
tration  cs  is  defined  at  the  stoichiometric  limit  of  L^Sis,  which  is 
equal  to  3.65  x  105  mol  m-3.  The  value  of  diffusion  constant  D  equal 
to  1 0-13  m2  s-1  is  taken  from  Ref.  [5].  It  can  differ  by  orders  of  mag¬ 
nitude  according  to  other  Refs.  [12,14,22,55].  The  Young’s  modulus 
of  copper  is  equal  to  128  GPa.  In  latter  simulations,  the  Si  thickness 
h  is  fixed  at  100  nm.  The  temperature  T  is  fixed  at  300  K. 

It  may  be  worth  noting  that  there  have  been  theoretical  and 
experimental  studies  suggesting  that  amorphous  Zintl  phase  com¬ 
pounds  occur  during  the  lithiation  and  delithiation  processes  of  Si 
regardless  of  its  initial  crystalline  or  amorphous  state  [56-60].  Thus, 
the  actual  variations  of  elastic  moduli  and  volume  can  be  more  com¬ 
plicated  than  what  are  depicted  in  Figs.  2  and  3.  This  is  also  why 


x/y  in  Li  Si 

X  y 


Fig.  3.  Variation  of  Li*Siy  alloy  volume  per  Si  atom  with  composition  ratio  x/y.  The 
dashed  line  indicates  a  linear  fitting  curve  to  the  data  points  (symbols). 


the  data  points  of  crystalline  Li— Si  alloys  were  connected  by  dotted 
lines  in  these  figures.  A  detailed  consideration  of  this  effect,  as  well 
as  of  the  aforementioned  possible  strong  composition  dependence 
of  Li-ion  mobility,  is  beyond  the  scope  of  the  present  work,  and  is 
left  to  a  future  study. 

3.  Finite  difference  method  for  stress-assisted  diffusion 


An  analytical  expression  of  stress  has  been  derived  above  for  a 
bilayer  beam,  with  a  Li-ion  concentration  field  as  a  loading  source. 
However,  the  diffusion  problem  is  nonlinear  and  can  only  be  solved 
numerically.  A  finite  difference  method  is  adopted  to  approxi¬ 
mate  the  transient  differential  equations  of  diffusion.  The  analytical 
expressions  of  stress  and  the  numerical  finite  difference  equations 
of  diffusion  are  combined  to  solve  the  coupled  diffusion  and  stress 
evolution  in  time.  The  numerical  scheme  is  summarized  below. 

The  Si  layer  is  discretized  in  /  equal  divisions  along  the  x  axis, 
as  shown  in  Fig.  4.  A  node  is  defined  at  the  middle  point  of  each 
division,  and  is  numbered  by  1, 2, . . .,  i, . . .,  /.  The  axis  of  time  is  dis¬ 
cretized  unevenly,  with  time  step  At  to  be  determined  by  current 
value  of  effective  diffusion  constant  D *  as  well  as  nodal  spacing  Ax 
(-/I//). 

The  mass  conservation  equation,  Eq.  (14),  is  approximated  by 

(18) 

where  l  indicates  the  Zth  time  step.  The  time  increment  At  is  deter¬ 
mined  by  a  A  x2/D*max,  where  D*max  is  the  maximum  nodal  value  of 
D  at  the  previous  time  step  /,  and  a  is  a  numerical  stability  control 
parameter,  set  to  be  a  =  1  in  later  simulations. 

The  flux  equation,  Eq.  (12),  is  approximated  by 
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where  I<  =  Mc(  1  -  c)(yE  -  (. EfC/E)crz)ky . 

The  boundary  conditions  at  i  =  1  /2  and  J  +  1/2  are  approxi¬ 
mated/given  by 

J (+1  =  A(c‘f  -  clf ),  (20) 

2 

=o.  (2i) 

l+2 

Note  that  in  the  special  case  that  the  Cu-layer  is  absent,  the 
problem  is  symmetric.  Only  one  half  of  the  Si  layer  needs  to  be 
considered,  and  the  same  boundary  conditions  as  above  (Eqs.  (20) 
and  (21 ))  can  be  applied  at  the  two  effective  ends. 
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The  solution  procedure  is  the  following.  Given  the  concentration 
field  at  the  Zth  time  step,  the  Si  modulus  and  its  derivatives  are 
calculated  according  to  Eq.  (17).  Then,  sz0  and  ky  are  obtained  by 
solving  Eqs.  (6)  and  (7),  which  leads  to  determination  of  the  stress 
field.  Then,  the  effective  diffusion  constant  D *  is  calculated  together 
with  the  bending-induced  flux  I<.  Finally  one  can  assemble  a  system 
of  algebraic  equations,  by  submitting  Eqs.  (19)-(21)  into  Eq.  (18), 
and  solve  for  nodal  concentration  c/+1(i  =  1,2, ...,  /)  at  a  current 
time  step  1+ 1.  The  process  is  repeated  to  obtain  the  diffusion  and 
stress  fields  in  time. 


4.  Numerical  results  and  discussion 

This  section  presents  the  simulations  based  on  the  method 
described  above  and  examines  the  effects  of  the  composition- 
dependent  modulus,  the  finite  concentration  (Eq.  (10))  and  the 
boundary  constraint  on  lithiation/delithiation  of  a  Si  and  a  Cu- 
coated  Si  nanorod.  Two  hundred  divisions  are  used  to  discretize  the 
anode  domain.  The  stability  control  parameter  a  is  set  to  be  one 
in  determining  each  time  increment  by  At  =  o'Ax2/D*max,  where 
D* max  is  the  maximum  effective  diffusion  constant  at  a  time  step. 
During  the  simulations,  ceq  is  initially  set  to  be  cs  for  the  lithiation. 
When  the  total  Li-ion  insertion  reaches  99%  of  the  total  saturation 
amount,  it  is  reset  to  be  zero.  Upon  that,  the  delithiation  process 
begins.  Other  than  this  switch,  the  solution  procedure  is  identical 
in  the  processes  of  lithiation  and  delithiation.  The  rate  constant  A 
in  the  boundary  condition  (Eq.  (15))  is  set  to  be  10-5  ms-1  for  all 
cases  examined.  The  numerical  convergence  was  achieved  with  the 
above  mesh;  doubling  the  mesh  density,  there  was  observed  little 
change  to  any  of  the  results  presented  in  the  figures  below. 


4.1.  Effects  of  composition-dependent  modulus 

At  first,  a  set  of  three  simulations  is  performed  in  the  absence  of 
a  Cu  layer.  Thus,  the  problem  is  symmetric,  and  the  Si  anode  stays 
straight  (i.e.,  no  bending).  In  this  case,  only  one  half  of  the  problem 
needs  to  be  considered,  as  discussed  before.  In  the  first  simulation, 
all  of  the  composition-dependent  modulus  E  and  its  derivatives 
E}C  and  ECc,  involved  in  the  above  diffusion  equations,  are  kept  in 
the  calculation.  In  the  second  simulation,  only  the  composition- 
dependent  modulus  E  is  kept,  but  all  of  its  derivatives  are  set  equal 
to  zero.  In  the  last  simulation,  the  modulus  E  is  set  to  be  a  con¬ 
stant  equal  to  the  modulus  of  Si.  Consequently,  all  of  its  derivatives 
are  equal  to  zero.  In  these  simulations,  the  total  Li-ion  amount  in 
the  Si  anode  is  recorded  over  time.  The  maximum  and  minimum 
stresses  in  the  anode  are  recorded  at  every  time  step.  Note  that 
the  maximum  and  minimum  stresses  shift  in  position  over  time. 
They  indicate  the  severity  of  stressing  that  the  anode  would  experi¬ 
ence  during  the  lithiation  and  delithiation  processes.  The  results  are 
plotted  in  Figs.  5  and  6,  respectively.  Fig.  5  shows  the  variation  of  Li- 
ion  uptake  over  time.  Fig.  6(a)  shows  the  variation  of  maximum  and 
minimum  stresses  over  the  lithiation  and  delithiation  processes. 
Fig.  6(b)  shows  the  variation  of  the  same  quantities  over  an  early 
short  period  of  time.  In  these  figures  and  all  other  figures  in  the  fol¬ 
lowing,  all  quantities  of  the  dimension  of  time  are  normalized  by 
characteristic  diffusion  time  tc  =  /i2/D.  All  quantities  of  the  dimen¬ 
sion  of  length  are  normalized  by  h.  All  quantities  of  the  dimension  of 
stress  are  normalized  by  Esi,  the  modulus  of  pure  silicon.  All  quanti¬ 
ties  of  the  dimension  of  concentration  are  normalized  by  saturation 
concentration  cs  at  the  stoichiometric  limit. 

From  Figs.  5  and  6,  it  can  be  seen  that  both  the  Li-ion  inser¬ 
tion  and  extraction  proceed  at  relatively  high  rate  initially  and 
then  gradually  slow  down.  It  results  in  high  stress  shortly  after  the 
beginning  of  lithiation  or  delithiation.  The  stress  decays  gradually 


Fig.  5.  Variation  of  Li-ion  uptake  (normalized  by  saturation  amount)  with  time.  Leg¬ 
end  “Full  E"  indicates  the  case  with  effects  of  composition-dependent  alloy  modulus 
fully  taken  into  account.  Legend  “Less  £  c”  indicates  the  case  with  alloy  modulus  vary¬ 
ing  but  its  derivatives  being  set  to  be  zero  in  the  diffusion  equations.  Legend  “Const 
£”  indicates  the  case  of  constant  alloy  modulus,  set  to  be  that  of  pure  Si. 

afterwards.  These  characteristic  behaviors  reflect  the  imposed 
boundary  condition. 

Compare  the  results  of  the  three  cases  of  fully  considered 
composition-dependent  modulus,  derivative-less  composition- 
dependent  modulus,  and  constant  modulus,  which  are  indicated  by 
“Full  E'\  “Less  E  c”  and  “Const  F\  respectively,  in  Fig.  5.  It  can  be  seen 
that  the  composition-dependent  modulus  plays  a  significant  role 
in  determining  the  diffusion  process  as  well  as  the  stress  field.  The 
effect  of  the  modulus  derivative  is  mild.  The  diffusion  curves  in  Fig.  5 
are  nearly  identical  in  those  two  cases  of  varying  modulus.  How¬ 
ever,  in  Fig.  6,  one  can  see  that  the  exclusion  of  modulus  derivatives 
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Fig.  6.  Variations  of  maximum  (solid)  and  minimum  (dotted)  stresses  with  time: 
(a)  full  lithiation  and  delithiation  process;  (b)  early  short  period.  Legend  “Full  £” 
indicates  the  case  with  effects  of  composition-dependent  alloy  modulus  fully  taken 
into  account.  Legend  “Less  £  c”  indicates  the  case  with  alloy  modulus  varying  but  its 
derivatives  being  set  to  be  zero  in  the  diffusion  equations.  Legend  “Const  £”  indicates 
the  case  of  constant  alloy  modulus  equal  to  that  of  pure  Si. 


B.  Yang  et  al.  /  Journal  of  Power  Sources  204  (2012)  168-176 


173 


Normalized  Time  tit 

c 


Fig.  7.  Variation  of  Li-ion  uptake  (normalized  by  saturation  amount)  with  time. 
Legends  “With  factor”  and  “Without  factor”  indicate  the  cases  with  and  without 
modification  factor  (1  -  c)  in  the  flux  equation,  respectively. 

in  the  calculation  makes  a  mild  difference  up  to  15%  in  the  pre¬ 
dicted  maximum  stresses.  This  is  in  contrast  to  the  estimate  of  1% 
difference  (in  potential)  made  by  Sethuraman  et  al.  [27]  We 
attribute  the  mild  difference  to  the  softening  composition- 
dependent  modulus  and  to  the  (elastic)  stress  magnitude  we 
allowed  to  reach  up  to  ~6%  of  Fsi  in  view  of  the  possible  high  yield 
strength  up  to  7-10 GPa  in  Si  [61,62]  compared  to  its  modulus  of 
1 70  GPa  and  the  measured  flow  stress  of  1 .5  GPa  in  Lii5Si4  by  Sethu¬ 
raman  et  al.  [25]  compared  to  the  measured  Young’s  modulus  of 
35  GPa  in  Li3Si  by  Sethuraman  et  al.  [26]  Note  that  the  significance 
of  the  terms  involving  F  c  is  determined  by  the  ratio  of  stress  to  the 
Young’s  modulus  (Eq.  (13)). 

4.2.  Effects  of  flux  equation  modification/finite  concentration 

Then,  a  simulation  is  run  where  the  classical  flux  equation  is 
applied  by  removing  factor  (1  -  c)  from  Eq.  (10).  All  other  consid¬ 
erations  are  kept  the  same  as  above.  The  same  results,  i.e.,  variations 
of  Li-ion  uptake  and  maximum  and  minimum  stresses  over  time, 
are  acquired  in  the  simulation,  and  are  plotted  in  Figs.  7  and  8. 

From  Figs.  7  and  8,  it  can  be  clearly  seen  that  this  modifica¬ 
tion  factor  (1  -c)  applied  to  the  classical  flux  equation  makes  a 
significant  difference  in  the  diffusion  process,  especially  in  the 
delithiation  process.  Although  the  predicted  stresses  are  about  the 
same  during  the  lithiation  process,  during  the  delithiation  process 
the  predicted  maximum  and  minimum  stresses  by  the  classical 
flux  equation  are  only  about  one  third  of  those  predicted  by  the 
modified  law.  This  is  expected  because  the  effective  diffusion  con¬ 
stant  is  much  lowered  at  high  concentrations  by  the  modification. 
With  the  same  rate  constant  of  insertion  and  extraction  at  the 
boundary,  the  diffusion  field  is  less  steep  and  thus  the  stresses 
are  lower  at  high  concentrations  within  the  classical  formulation. 
We  believe  that  the  present  formulation  with  the  modified  flux 
equation  is  more  physically  reasonable.  It  characterizes  the  early 
stages  of  lithiation  and  delithiation  in  a  unified,  symmetric  manner. 
It  predicts  the  early  processes  of  lithiation  and  delithiation  with 
about  the  same  magnitudes  of  stresses  and  about  the  same  dif¬ 
fusion  characteristics.  The  present  modification  symmetrizes  the 
description  of  the  lithiation  and  delithiation  processes.  If  these 
two  opposite  processes  are  different,  we  suggest  that  it  should  be 
caused  by  their  different  material  properties,  for  instance,  due  to 
different  mobility  of  Li-ion  migration  in  Si  and  in  Li44Si,  not  by 
the  law. 
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Fig.  8.  Variations  of  maximum  (solid)  and  minimum  (dotted)  stresses  in  Si  with 
time:  (a)  full  lithiation  and  delithiation  process;  (b)  early  short  period.  Legends  “With 
factor”  and  “Without  factor”  indicate  the  cases  with  and  without  modification  factor 
(1  -  c)  in  the  flux  equation,  respectively. 


4.3.  Effects  of  boundary  constraint-Cu  coating 

Finally,  a  set  of  simulations  is  performed  with  various  Cu-layer 
thicknesses  in  order  to  investigate  the  constraint  effect  of  a  Cu  layer 
bonded  to  one  side  of  the  Si  anode.  The  full  model  is  considered, 
including  the  effects  of  composition-dependent  modulus  and  its 
derivatives  and  the  flux  equation  modification.  Because  copper  is 
inactive  to  Li-ion  alloying,  the  Cu  layer  effectively  blocks  Li-ion 
from  flowing  through  the  coated  surface  area  of  the  Si  anode.  As 
discussed  before,  it  would  lead  to  the  anode  bending  and  hence 
asymmetric  diffusion  and  stress  fields.  The  variations  of  Li-ion 
uptake  and  maximum  and  minimum  stresses  in  the  anode  are 
recorded  over  time.  Also  the  variation  of  anode  bending  curvature 
is  recorded  over  time.  Since  the  characteristics  of  the  Li-ion  uptake 
variation  are  the  same  as  before  and  since  the  time  scales  involved 
can  be  viewed  from  the  other  profiles,  only  the  stress  and  bending 
curvature  variations  are  plotted  in  Figs.  9  and  10. 

Fig.  9(a)-(d)  shows  the  variations  of  maximum  and  minimum 
stresses  in  the  anode  over  time.  This  is  done  for  four  particu¬ 
lar  cases:  Cu-layer  thickness  equal  to  0  (trivial  thickness),  1  nm, 
lOOnm,  and  1  p>m.  Figs.  10(a)-(d)  shows  the  variations  of  anode 
bending  curvature  over  time  in  the  same  four  cases.  Recall  that  the 
anode  thickness  is  lOOnm.  The  Cu  modulus  of  128 GPa  is  roughly 
three  quarters  of  the  crystalline  Si  modulus,  and  roughly  three  and 
half  times  the  full  alloy  (L^Sis)  modulus.  From  these  figures,  it 
can  be  seen  that  the  time  scale  involved  in  the  lithiation  process 
is  very  sensitive  to  the  Cu  thickness  when  it  is  close  to  the  anode 
thickness.  The  time  to  finish  99%  of  the  lithiation  initially  increases 
with  Cu  thickness,  reaches  a  peak,  and  finally  decreases  with  further 
increase  of  Cu  thickness. 

When  the  Cu-layer  thickness  is  trivial  (Figs.  9a  and  10a),  the  Si 
anode  starts  to  bend  as  soon  as  the  lithiation  process  begins.  After 
a  short  period  of  time,  it  reaches  its  maximum  curvature  and  then 
gradually  recovers  from  bending.  The  Cu-layer  of  trivial  thickness 
blocks  the  Li-ion  flow  on  the  right  boundary  of  the  anode,  but  has 
no  mechanical  constraint  effect  on  the  anode  flexural  deformation. 
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Fig.  9.  Variations  of  maximum  (solid)  and  minimum  (dotted)  stresses  in  Si  with  time  at  various  Cu-layer  thicknesses:  (a)  0;  (b)  1  nm;  (c)  100  nm;  and  (d)  1  |xm. 


c 


Fig.  10.  Variations  of  anode  bending  curvature  with  time  at  various  Cu-layer  thicknesses:  (a)  0;  (b)  1  nm;  (c)  lOOnm;  and  (d)  1  |jim. 
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Fig.  11.  Variations  of  (a)  charging  and  (b)  discharging  times  at  various  percentages  of  completeness  (Li-ion  uptake  normalized  by  saturation  amount)  with  Cu-layer  thickness. 


The  resulting  nonuniform  concentration  field  causes  the  anode  to 
bend.  In  the  fully  charged  state,  the  concentration  field  returns  uni¬ 
formity,  and  the  anode  recovers  its  original  straight  configuration. 
The  maximum  and  minimum  stresses,  which  are  tensile  and  com¬ 
pressive,  respectively,  spike  at  the  early  time  but  quickly  decay 
along  with  the  lithiation  process.  When  discharging,  the  behavior 
is  characteristically  similar  to  the  above  case  of  charging,  except 
that  everything  switches  sign. 

When  the  Cu-layer  thickness  is  equal  to  1  nm  (Figs.  9b  and  1  Ob), 
which  is  one  hundredth  of  the  anode  thickness,  substantial  effects 
are  seen  in  the  stress  and  deformation  results.  The  anode  response 
is  about  the  same  as  above  in  the  early  period  of  time.  When 
the  lithiation  process  continues,  the  bending  curvature  and  the 
maximum  and  minimum  stresses  all  decrease.  However,  upon  full 
lithiation,  the  composite  structure  remains  bent  appreciably,  and 
the  resulting  stresses  in  the  anode  remain  high.  The  time  needed 
for  full  lithiation  is  more  than  doubled,  relative  to  the  above  case  of 
trivial  thickness.  In  contrast,  the  Cu  layer  of  thickness  equal  to  1  nm 
has  minor  effect  on  the  delithiation  process.  These  observations  can 
only  be  explained  by  invoking  that  the  1-nm  thick  Cu  layer  plays  a 
significant  role  in  both  of  the  diffusion  and  deformation  processes 
in  the  Si  anode  layer. 

When  the  Cu-layer  thickness  is  equal  to  lOOnm 
(Figs.  9c  and  10c),  the  diffusion  and  deformation  processes 
continue  to  change  relative  to  the  previous  cases.  The  time  needed 
to  fully  charge  the  anode  is  more  than  an  order  of  magnitude 
longer  than  that  in  the  first  case  of  trivial  Cu-layer  thickness.  The 
maximum  and  minimum  stresses  are  both  reduced,  indicating 
that  the  constraint  of  a  thicker  Cu-layer  holds  the  anode  in  more 
compressive  state  than  the  above  cases  of  nearly  free-standing 
anode  due  to  the  Li-ion  insertion.  The  bending  curvature  rapidly 
increases  initially  but  varies  very  slowly  afterwards,  indicating 
that  the  composite  structure  has  nearly  reached  a  lock-up  state 
of  deformation.  The  internal  stresses  only  vary  slowly  in  corre¬ 
spondence.  This  also  indicates  that  the  diffusion  process  evolves 
very  slowly  in  the  nearly  lock-up  state.  The  above  observations  are 
more  pronounced  when  the  Cu-layer  thickness  is  lOnm,  which  is 
also  examined  but  is  not  shown  here. 

When  the  Cu-layer  thickness  is  increased  to  1  p,m 
(Figs.  9d  and  lOd),  the  composite  structure  behaves  effectively  like 
a  thin  film  bonded  to  a  thick  substrate.  The  bending  curvature  is 
nearly  zero  during  the  entire  lithiation  and  delithiation  process. 
The  maximum  and  minimum  stresses  are  nearly  the  same  and 
are  both  compressive.  It  means  essentially  that  the  anode  is  fully 
constrained  by  the  Cu  layer.  Interestingly,  the  charging  time  is  on 
the  same  order  of  magnitude  as  in  the  case  of  trivial  constraint. 


In  order  to  better  understand  the  constraint  effect  at  var¬ 
ious  charging  and  discharging  stages,  the  time  needed  to 
charge/discharge  70%,  80%,  90%  and  99%  is  recorded.  This  is  done  for 
various  Cu-layer  thicknesses.  The  results  are  shown  in  Fig.  1  la  and 
b,  for  charging  and  discharging,  respectively.  It  can  be  seen  from 
Fig.  11a  that  between  1  nm  and  1  p>m  of  Cu-layer  thickness,  the 
charging  time  increases  substantially,  especially  for  high  charging 
percentages.  The  peak  can  be  a  few  orders  of  magnitude  longer  than 
that  of  either  trivial  or  very  thick  Cu-layer  support.  In  contrast,  the 
Cu-layer  constraint  effect  has  relatively  little  effect  on  the  discharg¬ 
ing  process  in  the  entire  range  of  examined  Cu-layer  thicknesses. 
By  examining  Eq.  (12)  and  Fig.  11a  and  b,  it  may  be  remarked  that 
the  term  independent  of  concentration  gradient  but  proportional 
to  bending  curvature  in  Eq.  (12)  is  always  opposing  to  the  Cu  layer. 
In  the  lithiation  process,  it  can  lead  to  a  lock-up  state,  effectively 
blocking  the  anode  from  full  charging.  In  contrast,  in  the  reversed 
delithiation  process,  this  term  assists  the  Li-ion  diffusion  out  of  the 
anode,  which  is  clearly  seen  in  Fig.  lib. 

5.  Conclusions 

A  fully  coupled  chemoelastic  analysis  has  been  carried  out  for 
lithiation  and  delithiation  of  a  bilayer  Cu-coated  Si  anode.  The  clas¬ 
sical  Larche-Cahn  chemical  potential  [42]  was  adopted  to  describe 
the  Li— Si  alloy  system.  The  term  of  chemical  potential  due  to  elas¬ 
tic  modulus  variation  with  composition  is  kept  in  order  to  evaluate 
its  effect  on  Li-ion  diffusion.  A  nonlinear  flux  equation  is  intro¬ 
duced  that  symmetrizes  the  law  for  Li-ion  diffusion  at  dilute  and 
near-saturation  concentrations.  Meanwhile,  the  deformation  and 
stress  chemically  induced  by  Li-ion  insertion  in  the  bilayer  com¬ 
posite  structure  are  modeled  within  the  classical  beam  theory  and 
solved  analytically.  The  stress-assisted  diffusion  problem  is  non¬ 
linear  and  is  solved  numerically  by  applying  a  finite  difference 
method. 

It  has  been  shown  that  the  effects  of  composition-dependent 
Li— Si  alloy  modulus  are  significant  on  the  stress  field  but  less  sig¬ 
nificant  on  the  diffusion  process.  The  second  derivative  of  modulus 
with  respect  to  composition  can  be  safely  omitted  from  the  dif¬ 
fusion  equation.  In  contrast,  the  first  derivative  can  make  a  mild 
difference  up  to  1 5%  in  predicted  maximum  stresses  with  and  with¬ 
out  it.  Also,  it  has  been  shown  that  the  flux  equation  modification 
can  make  a  significant  difference,  especially,  during  the  delithiation 
process.  The  modified  flux  equation  is  considered  as  more  reason¬ 
able  when  applying  to  the  whole  range  of  concentration  from  dilute 
to  near  saturation.  It  describes  the  Li-ion  diffusion  at  dilute  and 
near-saturation  concentrations  in  a  unified,  symmetric  manner. 
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Finally,  it  has  been  shown  that  the  Cu  layer  on  one  side  of  the  Si 
anode  can  play  a  profound  role,  not  only  in  blocking  surface  Li- 
ion  flux,  but  also  in  substantiating  a  term  of  bending  curvature 
in  the  flux  equation  opposing  Li-ion  insertion  into  the  anode.  This 
term  due  to  bending  stress  gradient  can  effectively  hold  the  anode 
from  full  lithiation  for  a  range  of  Cu-layer  thickness,  but  can  also 
somewhat  accelerate  the  delithiation  process. 
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